Further experimental results on modelling and algebraic control of a delayed looped heating-cooling process under uncertainties

The aim of this research is to revise and substantially extend experimental modelling and control of a looped heating-cooling laboratory process with long input-output and internal delays under uncertainties. This research follows and extends the authors' recent results. As several significant improvements regarding robust modelling and control have been reached, the obtained results are provided with a link and comparison to the previous findings. First, an infinite-dimensional model based on mathematical-physical heat and mass transfer principles is developed. All important heat-fluid transport and control-signal delays are considered when assembling the model structure and relations of quantities. Model parameter values optimization based on the measurement data follows. When determining static model parameter values, all variations in steady-state measured data are taken into account simultaneously, which enhances previously obtained models. Values of dynamic model parameters and delays are further obtained by least mean square optimization. This innovative model is compared to two recently developed process models and to the best-fit model that ignores the measured variations. Controller structures are designed using algebraic tools for all four models. The designed controllers are robust in the sense of robust stability and performance. Both concepts are rigorously assessed, and the obtained conditions serve for controller parameter tuning. Two different control systems are assumed: the standard closed-loop feedback loop and the two-feedback-controllers control system. Numerous experimental measurements for nominal conditions and selected perturbations are performed. Obtained results are further analyzed via several criteria on manipulated input and controlled temperature. The designed controllers are compared to the Smith predictor structure that is well-established for time-delay systems control. An essential drawback of the predictor regarding disturbance rejection is highlighted.


Introduction
Heating-cooling loops equipped with a heating power source and a heat sink connected by pipes appear typically in engines used in fuel-based power plants [1,2], and in cooling systems of an automotive combustion engine [3,4]. Their complete structures adopt a wide range of possible topologies [5]. The heat source and sink components can generally be viewed as heat-exchangers (HXs). Numerous academic research, experimental and industrial applications of processes with HXs and their networks have been of a heating-stations pipe network system to reduce energy consumption and carbon emission was proposed in Ref. [15]. Pipe network and building inertia thermal delays were included in the design.
In a recent work [38], an adaptive energy optimization control technique with a disturbance estimator, an energy consumption optimization block, a state predictor, and a state tracking controller for an advanced engine thermal management system has been proposed. The authors have pointed out that some control strategies [39,40] for looped coolant systems with a long dead time cause a large overshoot in temperature control. However, the internal delays were not taken into consideration. Moreover, the presented experiment has resulted in a non-zero reference-tracking steady-state error.
The robust control framework attempts to design and tune a control system sufficiently insensitively to process perturbations, model uncertainties, and signal disturbances. Especially, control system stability and a desired performance level should be satisfied. Besides advanced and computationally demanding robust control strategies applied to looped heating-cooling systems without considering delays in the design explicitly [11,41,42], other methods and results incorporating delayed models have been published. A robust temperature controller of a fluid-fluid HX process with uncertainty estimation was proposed in Ref. [43]. Therein, a delay induced by the actuator dynamics was considered in simulation tests. Santos et al. [14] designed a model predictive controller based on a quadratic program solution for a nonlinear solar collector model, where a robust input-output delay compensation scheme was proposed. A combination of the PI controller and the Clegg integrator compensator in the two-step robust temperature control under uncertainties and dominant time delay effect on the reset action was investigated in Ref. [44]. The controlled experimental HX system used in the food industry was modelled by a set of first-order plus time delay submodels. Using the same dead-time submodel, six robust controllers in real-time controlling a laboratory HX with nonlinear and asymmetric dynamics and with process gain, time-constant, and time-delay uncertainties were compared in Ref. [45]. The second-order plus time delay model of an HX process was used for robust controller design, ensuring stable control performance with model uncertainties in Ref. [46]. The controllers were tuned by the celebrated Ziegler-Nichols rules and by H-infinity synthesis. The third-order plus time delay model of a shell-and-tube HX was applied to derive, verify, and compare several robust controllers in Ref. [47]. The controllers were benchmarked using the integral absolute error (IAE) criterion. A model of a plate HX with a nonlinear static characteristic and time delay was utilized when performing a robust fuzzy evolving cloud-based control design [48]. Gupta et al. [16] considered different uncertainties (including that in the input-output delay value) in a linear model of a liquid-liquid HX process when deriving parameters of the conventional PID controller using the H ∞ robustness metric.
It is worth noting that input-output delay models only (i.e., neither DDEs nor DPDEs) were considered when performing the abovegiven research. Results on advanced control of looped HX processes incorporating delayed infinite-dimensional models [5,[28][29][30], however, have not utilized robustness methods and tools in the design explicitly. Moreover, they are mostly mathematically and computationally demanding. Recall that infinite dimensionality is natural for HXs and their loops; the former is due to an infinite-dimensional solution of PDEs, the latter because of the delay effect of long piping.
In [39], a control design method based on using algebraic tools [49] in the input-output space satisfying robust stability and performance was proposed when controlling the fluid temperature in a laboratory looped heating-cooling process. A Two-Feedback-Controllers (TFC) control system structure was used therein. A nonlinear process model derived based on DDEs was further linearized in the vicinity of a steady-state operating point for control design aims. Note that the results have been slightly refined in Ref. [50] later. The model and its parameters were obtained using a three-step procedure [26]. The second step estimated static parameters for each single process submodel separately and the results of each substep were used to determine the parameters of the following submodel. Another approach to identifying the model parameters based on the relay-feedback experiment was published in Ref. [27]. While the computational and experimental burden was less than in the previous method, worse results were obtained. However, the authors hypothesized that despite a low accuracy, the eventual model could be sufficient for control purposes.
Hence, the presented research has been motivated by gaps in the modelling and control design in the above-mentioned results and research questions raised therein. Moreover, the research is driven by a constant effort to develop advanced control techniques to reduce energy dependence while maintaining sufficient user comfort.
The motivation can be summarized as follows: a) Using the mathematical-physical modelling via DDEs, the complete set of static parameters can be identified at once employing the least-square technique, i.e., not for each submodel separately. The obtained model should better cover and interpolate all the measurement variations and uncertainties than the original model [26], which adopted a part-by-part procedure. b) Another two process model parameter sets are considered for the control design. First, it is a model that best estimates the nominal (unperturbed) process response. Such a model might perform worse robustness despite its better nominal accuracy. Second, the best model obtained by the relay-feedback test [27] is assumed. The research question is whether such a model can yield sufficient control responses. c) The use of the standard simple One-Degree-of-Freedom (1DoF) control system can be compared to the TFC structure. The goal is to track the step-wise and linear-wise reference temperature value and asymptotically reject the constant load disturbance using algebraic tools. d) In Ref. [39], the robust performance (that also includes robust stability) was not evaluated in detail. Therefore, the reader should be provided with detailed evaluations of all four models. e) Both the control systems combined with all four models are benchmarked when real laboratory measurements on the appliance.
The nominal case and three perturbed cases are considered. f) Control target of engine cooling loops is to reach desired temperature levels while it is to reduce heat consumption in HX networks [5]. Hence, we consider both criteria when evaluating the experimental results. g) As the process model structure includes input-output delay, the thing to ask is to use the well-established Smith predictor dead-time compensator [51] as an alternative to 1DoF and TFC. Hence, the Smith predictor controller and its properties are derived.
The rest of the paper is organized as follows. Section 2 concisely introduces the experimental laboratory appliance and its mathematical model via DDEs and auxiliary algebraic relations. Optimized robust model parameters identification is presented in Section 3. In addition, the unperturbed best-fit model, the original model [26], and the relay-based model [27] are given to the reader in the section as well. Controller structures via algebraic means for 1DoF and TFC control systems are derived in Section 4. Section 5 provides robust stability and performance conditions and their application to the derived models and control structures, giving rise to controller parameter settings. Real laboratory experiment results are displayed and evaluated by several performance measures in Section 6. Section 7 concisely links to the Smith dead-time compensator and highlights its drawback. Then the paper is concluded, and possible future research is sketched.
Note that only necessary information is placed in the main text body to be concise. Previously published data and extensive mathematical derivations and proofs are given in appendices, or the reader is referred to the literature. The basic used notation is summarized in a table.

Laboratory appliance and process model structure
A concise description of the looped heating-cooling process and the used laboratory test bed follows. A process model based on DDEs arising from heat balances and some auxiliary relations is also given to the reader [26,39].

Looped heating-cooling process description
The looped heating-cooling process is sketched in Fig. 1 [39]. The process works as follows [26]: The heating fluid (distilled water here) in the loop is driven by a centrifugal pump controlled by the input voltage u P (t). It yields a change in the fluid flow rate ṁ(t) in the piping. The fluid flows through an instantaneous water heater, the input power of which is P H (t). The heater outlet temperature of the fluid is ϑ HO (t). The fluid then goes via a long well-insulated pipeline into a solid-liquid plate-and-fin HX (i.e., a radiator) that serves as a heat sink (let us call it a "cooler" for simplicity). The cooler fluid inlet temperature is measured as ϑ CI (t). The rate of heat flow depends on the ambient temperature ϑ a (t). The cooler fluid outlet temperature ϑ CO (t) is affected mainly by the input voltage u C (t) of a fan connected to the cooler. The voltage value can be controlled continuously, whereas the second fan can only switch on or off. Changes in the heat fluid volume are compensated by a small expanse tank placed close to the cooler. Finally, the fluid flows into the pump, which closes the loop.

Laboratory appliance appearance and equipment
The front and back sides of the laboratory appliance realizing the looped heating-cooling process are displayed in Fig. 2. Its technical description follows as per the positions indicated in the figure. Note that the reader is referred to Ref. [26] for further details about electronics inside the laboratory model, wiring, connection to a PC, and HW and SW equipment on the PC side.

DDEs-based process model structure
The modelling procedure based on the deductive principle has three steps [26]: Step 1. Every distinct functional part of the process is modelled separately using DDEs based on heat transfer balances and via auxiliary algebraic equations expressing static relations of particular quantities. The obtained submodels are then linked via their mutual physical quantities due to fluid flow in the piping, which gives rise to delays. The crucial idea here is that all significant delays are considered.
Step 2. Heat transmission coefficients, variables in algebraic relations, and the mass flow rate are identified by steady-state inputoutput measurements, i.e., all derivatives in DDEs are assumed to be zero.
Step 3. Masses and delays are determined via measured dynamic responses by minimizing integral criteria in the time domain.
The following equations in this subsection summarize Step 1, whereas, Steps 2 and 3 are solved anew in Section 3 to improve the model robustness. Note that the following assumptions are taken into account when modelling: water inside the pipelines and functional parts of the process is incompressible; the fluid flow ṁ(t) is constant in any place of the process within a particular time instant; the specific heat capacity c (≈ 4175 J kg − 1 K − 1 ) of the fluid is constant within a range of operating temperatures; the specific heat capacity of the piping wall material (copper) is neglected; the dynamics of sensors and actuators is omitted (see also Remark 1).
The heater dynamics is given by the imbalance between the inlet and outlet fluid heat, the input heat power from the solid surface of the heater, and the waste heat leaking through the heater shroud to the ambient environment. The particular DDE reads where M H means the water mass inside the heater and k H (t) stands for the heat transmission coefficient, the value of which is assumed Laboratory appliance with the looped heating-cooling process [26].
to be dependent on P H (t) and ṁ(t) as where h i , i = 0, 5 are real-valued constants.
The duration of the heat fluid flow through the heater is τ H . Hence, the delay value 0.5τ H in (1) expresses that the heat power acts in the middle of the heater. The long pipeline between the heater and the cooler can be modelled based on the heat balance and asynchronous effect of fluid temperatures on the pipeline inlet and outlet as where M P means the water mass inside the long pipeline, k P denotes the heat transmission coefficient (that is considered being a constant because of its low value), and τ HC means the time during which the fluid goes from the heater outlet to the cooler inlet. Note that the arithmetical mean temperature on the right-hand side of (3) agrees with the case that the heat loss is constant along the pipeline at a time instant [26]. The cooler dynamics is modelled by a DDE analogously to (1) and (3) as where M C means the water mass inside the cooler plates and tubes, τ C denotes the fluid-flow transport delay between the cooler inlet and outlet, and k C (t) stands for the heat transmission coefficient, the value of which is assumed to be dependent on delayed control voltage u C (t).
where τ FC is the fan reaction delay taking into account latencies of electronic and mechanical fan parts, and c i , i = 0, 2 are real-valued constants. Note that the on/off fan (that influences k C (t)) is permanently switched on. The mass flow rate is related to u P (t) as follows in the model where π i , i = 0, 2 are real-valued constants. The fluid-flow delay τ CH between the cooler outlet and the heater inlet (going through the pump) is modelled via a static relation only because of a short distance that means that heat loss is omitted. shapes of unit step responses on ΔP H (t) (close to steady states s ϑ ⋅ ) are taken as representatives in the figure.
The following notation is used in Fig. 1: Remark 1 Sensor delays are not considered in the model for several reasons, most of which yield from the fact that model delays are estimated based on the measurement, not derived analytically. First, as can be deduced from Fig. 1, almost all delays (except for τ H and τ FC ) are relative to time instants when significant shape changes of the characteristics appear. It is sufficient to detect the changes, not the temperature values themselves, for the delay values identification. Hence, these changes do not depend on sensors' time constants (≈ 8 s) that are much less than the overall internal delay (≈ 140 s) due to the loop closeness. In addition, possible sensor latencies have only a minor effect on the overall dynamics since they act in inputoutput relations only [27]. Last but not least, correct values of dead times are not crucial when experimental controlling the process because both delay-evaluation characteristics and controller inputs adopt the same measurement data.

Model parameters identification
The model parameters identification (see Steps 2 and 3 introduced in Subsection 2.3) is thoroughly revised and reformulated in this section, compared to the original result given in Ref. [26]. Recall that parameter values of every single functional part of the process were estimated separately therein. Now, let us perform the innovative concept that considers model (1), (2), (3), (4), (5), (6), and (7) at once.
Consider the following ranges of inputs The range of u P has been chosen relatively narrow for two reasons. First, the proposed control strategy is valid for constant delays; therefore, the flow rate induced by the pump input voltage should not significantly vary, and it ought to stay near the operating point. Second, u P does not pose to be the manipulated input. A relatively wide range of u C can serve to verify the robustness properties of the proposed control system by performing perturbations. The range of P H values is also wide as it serves as the manipulated input for controlling particular fluid temperatures. Besides, these heating power values are not symmetric with respect to the operating point but are shifted up. This feature is because of the assumption that heating rather than cooling is required for the process when controlling fluid temperature. It is also worth noting that the static characteristics for slightly wider ranges of u P , u C , and P H can be estimated using linear extrapolation.
Steady values of fluid temperatures have been measured for ranges (9), as summarized in Appendix A [26]. Note that the values are normalized to the nominal constant ambient temperature of ϑ a = 24 • C. Substitute (2), (5), and (6) into (1), (3), and (4). Set zero derivatives on the left-hand sides of DDEs and let all variables be at their steady states (i.e., the delays vanish as well). Then, the set of obtained three nonlinear algebraic equations can be expressed in a condensed form where 0 = (0, 0, …, 0) T , v T = ( s u P ,s u C , s P H , s ϑ HO , s ϑ CI , s ϑ a ), and p T = (h 0 ,...,h 5 ,π 0 ,π 1 ,π 2 ,c 0 ,c 1 ,c 2 ,k P ). Define by J(v, p) the 3 × 13 Jacobian of f(v, p) with respect to p, i.e., each expression on the left-hand side of (10) is subject to partial derivatives according to each element in p.
Each line of measured data in Table A1 constitutes one real-valued vector 34 and 34 . Hence, f(p) and J(p) have dimensions of 102 × 1 and 102 × 13, respectively.
Then, the Levenberg-Marquardt method [52] enables finding the solution of f(p) = 0 iteratively via (11) where k means the iterative step. The method represents a stochastic approximation technique that minimizes the H 2 norm of f(p), i.e., it solves nonlinear least-square problem (10).
Since the solution p * = lim k→∞ k p significantly depends on the initial estimate 0 p and the evolution of the damping factor γ > 0, it might not be unique. Table 1 displays the two best solutions obtained numerically, including the original result [26] for comparison.
We call the mathematical model (1) to (7) with parameter values in the rightmost row of Table 1 the "Original model" hereinafter. As can be seen from the table, the approximation error measured by the H 2 norm has been reduced compared to the Original model. Another advantage is that the results cover a broad spectrum of process perturbations and measurement uncertainties, which is beneficial for the further robust control design. On the other hand, it may lead to a worse estimation of the nominal case.
Nevertheless, Result 1 in Table 1 is unsuitable for robust control design since it enables only a reduced range of pump voltage input (i.e., a range of fluid flow values). Namely, it is clear from (6) that whenever u P (t) < − π 1 , the flow rate ṁ(t) becomes complex-valued.
That is, Result 1 admits u P (t) ≥ 3.7155 V, which is rather limiting (regardless, it complies with ranges (9)). Therefore, we have selected Result 2 in Table 1 for further identification and control in this research.

Dynamic parameters estimation
The so-called "dynamic" parameters influence the transient part of a time-domain response (see Fig. 3). These are masses M H , M P , M C and all delays τ H ,τ HC ,τ C ,τ CH ,τ FC . The error minimization between the measured and modelled responses can estimate their values.
However, as the static parameters determined in the preceding subsection lead to incorrect static gain estimation (due to perturbations), this gain can be easily adjusted by an additional gain, so that the steady-state values coincide. Nevertheless, such an adjustment may harm other decisive parts of the responses indicated in Fig. 3. Therefore, the dynamic parameter values have eventually been found for Result 2 of Table 1 based on the following intuitive assumption: Assumption 1. Consider a model that matches the nominal responses perfectly. Then, its dynamic parameter values must be optimal for the model under perturbations.
To rephrase Assumption 1, the complete model that fits the measured nominal measured data in the vicinity of some operating point optimally is found first. Let us denote this model simply as "Best-fit". Then, the model partially found in Subsection 3.1 adopts the dynamic parameters of the best-fit one.
As the aim is to use P H (t) as the manipulated input for the control tasks, only responses (of temperatures) to the step change ΔP H (t) are assumed when finding the best-fit model parameters. Let the objective be the IAE between the measured (ϑ ⋅ (t)) and modelled (ϑ ⋅,m (t)) responses where Δt is the sampling time. Hence, the cost function of the optimization problem for the best-fit model can be defined as follows: Let us select the following operating point for the nominal data The solution of (13) for t 0 = 0 s, t 1 = 2000 s, Δt = 1 s using the well-established Nelder-Mead flexible-simplex method [53] yields the dynamic parameters provided in Table 2. The table also contains the values of the original model [26].
The value of τ FC have been found graphically based on the response to Δ s u C (t). Note that other delays can alternatively be found based on Fig. 1. That is, the four particular time intervals can also be identified from the measured data, giving rise to τ H ,τ HC , τ C ,τ CH .
Let us denote by "Model 1" the novel model having static parameters introduced as "Result 2" in Table 1 and dynamic parameters provided as "Best-fit" in Table 2.
The knowledge of nominal delay values and particular modelled relations u P (t) ↦ṁ(t) enables estimating the operating ranges (or possible perturbations) of the delays. Relations (15) express that delays are in inverse proportion to ṁ(t).
The obtained ranges are summarized in Table 3.

Model linearization and results comparison
As the aim of this research is to apply robust control principles for linear systems, nonlinear model (1)- (7) has to be linearized. Appendix 2 follows a simple linearization procedure in the neighbourhood of the operating point (14) introduced in Ref. [26].
A comparison of unit step responses for input and output variables is displayed in Fig. 4(a-c). Original model, Model 1, and the Bestfit model are included for the benchmark. Fig. 5 includes those unit step responses for the unified static gain. As mentioned above, the modelled responses can simply be adjusted using an auxiliary gain in practice.
The figures prove the Best-fit model superiority when matching the nominal step responses. Plots for the Original model and Model 1 in Fig. 4(a-c) evince significant errors when responding to Δ s u P (t) and Δ s u C (t). The main reason is that both models attempt to cover the perturbations and uncertainties, and the nominal case represents only a limited subset of the measured data. It is also partially due to the dynamic model parameters have been found based on the response to Δ s P H (t); however, transient parts of responses are estimated quite well, as clear from Fig. 5(a-c). Whereas the Original model works better in cases (b) and (c), Model 1 matches the measured responses better in case (a). It must, however, be highlighted again that the best matching of the nominal case is not the primary goal of these models.
In [27], various models of the looped heating-cooling process in question were received by a relay-feedback experiment. These relay-based models have been computed for relation P H (t) ↦ ϑ CO (t) only. It can be deduced from (B.1) that the relation is given by the DDE where τ 0 = 0.5τ H , τ a = τ H + τ HC + τ C + τ CH , τ b = τ HC + τ C . Four particular models with the best integral measures in the time and frequency domains are provided in Table 4. Coefficients of (16) for the Original model, Model 1, and the Best-fit model are also included in the table for completeness. Note that this data will be used in the next sections. Surprisingly, although Model 1 adopts the dynamic parameters of the Best-fit model, the coefficients of the latter one are closer to the Original model. Let the best model (measured by the IAE) in the first column (i.e., Relay 1) of Table 4 be the nominal "Relay model" for control purposes. Notice that non-delay coefficients differ from other models significantly, and the internal delay τ a is far from the physical one of the process. The reason for including the relay-based models is the following. This research question was raised in Ref. [27]: Can the process be sufficiently controlled based on the (relatively inaccurate) relay-based model? Hence, a partial goal of this research is to benchmark the model usability in the robust control design against other process models. Unit step responses Δ s P H (t) ↦ Δ s ϑ CO (t) for the Relay model with/without the static gain adjustment are displayed in Fig. 6. Note that the IAEs (after the unifying of the static gains, i.e., the steady-state temperatures coincide) for the Original model, Model 1, Best-fit model, and Relay model are, respectively, 0.114, 0.098, 0.053, and 0.984.
A frequency-domain model comparison is made via Nyquist plots; see Fig. 7(a-i) for unadjusted static gain and Fig. 8(a-i) for adjusted ones. Note that the actual measured process data was obtained simply by inserting the sinus signal into the particular input (at the operation point) and observing the amplitude and phase shift, not via the Fourier series analysis. In the figures, G ij (s) means the entry of the transfer function matrix G(s) in (B.1).
Unadjusted nominal frequency responses for the Best-fit model naturally provide the best matching among all the models. Although the shapes of all plots are similar, their scales differ significantly, mainly due to process perturbations that have to be covered. Contrariwise, Nyquist plots with adjusted static gains are similar in shape and size in many cases. In Fig. 8(a-i), the Original model and Model 1 estimated the measured data well (besides the Best-fit model does).

Controller structure design by algebraic tools
In this section, the reader is provided with a concise description of the used control system structures, namely, 1DoF and TFC. Parameterized control laws are then derived for each structure using algebraic tools introduced in Appendix C.
The goal is to control ϑ CO (t) using the manipulated input P H (t) based on submodel (16). Other external inputs, i.e., u P (t), u C (t), ϑ a (t) are considered being model uncertainties. The transfer function of relation (16) reads

1DoF and TFC control systems
The 1DoF control system agrees with the simple negative control feedback loop; see Fig. 9. In the figure, r(t), e(t), u(t), d(t), y(t) denote the reference, control error, control action (i.e., the computed manipulated input), load disturbance, and the output variables, respectively. The controlled system and controller transfer functions are G(s),C R (s), respectively.
The TFC scheme displayed in Fig. 10 includes an additional inner feedback loop that may help to stabilize the control system, adjust its dynamics or attenuate the disturbance effect by an additional degree of freedom. The transfer function of the controller acting inside the inner loop is denoted as C Q (s).

Algebraic representation of signals and transfer functions in the control systems
Algebraic control design in the ring of quasi-polynomial meromorphic functions (R QM ) is adopted [49,50]. The reader is referred to Appendix C for the ring definition (Definition C.4). The controller design idea is based on the fractional representation of all control system transfer functions and signals using R QM . One can write The following three goals to be satisfied are as follows: a) The control system is stable, i.e., all transfer functions are in R QM . b) The output asymptotically approaches the reference, i.e., lim t→∞ y(t) = lim t→∞ r(t) or lim t→∞ e(t) = 0 for d(t) = 0, r(t) ∕ = 0. c) The load disturbance is asymptotically attenuated, i.e., lim t→∞ y(t) = 0 for d(t) ∕ = 0, r(t) = 0.

Controller structure design for 1DoF
The following three lemmas can be proven [49].
Note that a particular solution pair R p (s), P p (s) can further be parameterized as for P(s) ∕ = 0, Z(s) ∈ R MS .

Lemma 2. The reference is asymptotically tracked if and only if (21) holds
Now, consider controlled process model (17). A controller that satisfies all three above-given conditions for a linear-wise r(t) and a step-wise d(t) (that is a natural demand in practice) has the transfer function C R (s) = m 0 m 1 a(s)(r 1 s + r 0 ) p num (s) where p num (s) = p num,4 s 4 + p num,3 s 3 + p num,2 s 2 + p num,1 (s)s + p num,0 (s), p num, and m 0 , m 1 > 0 are selectable controller parameters. A detailed controller derivation is given in Appendix D.

Controller structure design for TFC
For the TFC structure, the following lemmas hold [39,50].
where V(s) can be written by (26) A particular solution V p (s), P p (s) can further be parameterized as for P(s) ∕ = 0, Z(s) ∈ R MS .

Lemma 5. The reference is asymptotically tracked if and only if
Consider a linear-wise r(t) and a step-wise d(t) again. Possible controllers satisfying control feedback stability, asymptotic reference tracking, and load disturbance rejection are governed by the transfer function C Q (s) = m 3 0 a(s)λv 1 s 2 p num (s)(s + m 1 ) , where p num (s) = p num,4 s 4 + p num,3 s 3 + p num,2 s 2 + p num,1 (s)s + p num,0 (s), p num, and m 0 , m 1 > 0, λ ∈ (0, 1] are selectable controller parameters. Clearly, the TFC control system has one more tunable parameter than the 1DoF system (in our case). Moreover, reference tracking and disturbance rejection tasks can be partially solved independently. A detailed derivation of the controllers can be found in Appendix E.

Robust stability and performance
Based on models obtained in Sections 2 and 3, parameters of controllers (23)- (24) and (30)- (31) are set to meet robustness requirements. Namely, robust stability and robust performance are considered [54]. Roughly speaking, when these requirements are met, the control feedback system remains stable and provides a satisfactory response under process perturbations and model uncertainties. The perturbations are due to all variations of uncontrolled inputs, see (9), and the fluctuation of the ambient temperature ϑ a ∈ [18,28] • C. Recall that the nominal setting is given by the operating point (14).
The nominal model transfer function G m,0 (s) have structure (17) with parameters given in Table 4. The set of all perturbed transfer functions G m (s) satisfies unstructured multiplicative uncertainties (32) where ‖Δ(s)‖ ∞ ≤ 1 is a stable bounded function. W M (s) expresses a fixed stable weight function of the uncertainty frequency distribution and is searched so that inequality is satisfied without an excessive conservativeness. Definition 1. The control system is robustly stable if it is stable for G m,0 (s) and remains stable also for all G m (s).

Definition 2. The control system satisfies the nominal performance if it holds that
where W P (s) is the sensitivity weight function and S 0 (s) stands for the nominal sensitivity function that agrees with the transfer function between r(t) and e(t) for the nominal model G m,0 (s).
To rephrase Definition 2, inequality (34) expresses that the nominal control system performs "sufficiently well", as per the selected bound 1/W P (s).
where S(s) means the sensitivity function for the whole set of G m (s).
Particular conditions for 1DoF and TFC control systems and their application to the heating-cooling process models and their derived controllers follow.

Robustness design for 1DoF
The following lemmas hold for the 1DoF control system under the assumption that G m,0 (s) and G m (s) have the equal number of poles s i with Res i ≥ 0 [54].

Lemma 7. The control system is robustly stable if and only if
where T 0 (s) = 1 − S 0 (s) is the nominal complementary sensitivity function that agrees with the transfer function between r(t) and y(t) for the nominal model G m,0 (s). For 1DoF with controller (23)- (24), it holds that with the nominal parameters in Table 4.

Lemma 8. The control system satisfies inequality (35) and robust stability condition (36) if and only if
Note that S 0 is given by (39) Hence, the estimation of W M (s) followed by the selection of W P (s) and testing inequality (38) for the Original model, Model 1, the Best-fit model, and the Relay model is presented below.  Table 3, and the ambient temperature range ϑ a ∈ [18,28] • C are considered (color lines). The eventual upper bounds (40) are displayed (thick black line). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Fig. 11(a-d) displays the left-hand sides of (33) upper-bounded by |W M (jω)| that are constructed by virtue of the procedure described in [39] (see Section 4 therein). Perturbations (9), the corresponding delay variations in Table 3, and the ranging of the ambient temperature ϑ a ∈ [18,28] • C are taken to calculate the set of all G m (s). For the Relay model, it is impossible to apply the above-given perturbations. Models "Relay2", "Relay3", and "Relay4" in Table 4 are taken as modelled perturbations in this case.
The following idea is used in the setting of W P (s) in (34). Let s i = − a be the dominant pole of G m,0 (s). It is generally suggested to set it as the pole of T 0 (s). I.e., assume the approximation T 0 (s) ≈ a/(s +a) that agrees with S 0 (s) ≈ s/(s + a). That is, one can consider the upper bound in (34) as W P (s) ≈ (s + a)/s. However, numerical tests have shown that it is impossible to meet the condition; therefore, a weaker requirement is eventually taken, with a sufficient margin on (34). Hence, it is set m 0 = m 1 = a for the nominal model controller (23)-(24) and k W is found so that ‖S 0 (jω)‖ ∞ = x‖W P (jω)‖ − 1 ∞ where x ≈ 0.6 is the selected margin (conservativeness).  Fig. 12(a-d).
The Original model returns the minimum value of 0.930 for m 0 = 0.003, m 1 = 0.018 and the Relay model that of 0.615 at m 0 = 0.002, m 1 = 0.011, which represents satisfactory results. However, Model 1 returns the minimum value slightly above the performance border. The Best-fit model gives values much higher than 1. Therefore, two additional tests are performed in Fig. 13(a and b) In none of the two tests, a significant improvement has been obtained. Model 1 gives 0.987 for m 0 = 0.003, m 1 = 0.01 that is too close to the performance border yet satisfactory. The Best-fit model remains robustly unstable; hence, the robust performance cannot be met (the minimum is 2.198 at m 0 = 0.006, m 0 = 0.2). Besides, higher values of m 0 , m 1 mean faster control responses that usually imply excessive overshoots. The eventually selected controller parameters are summarized in Table 6.

Robustness design for TFC
The following lemmas hold for the TFC robust control design under the assumption that G m,0 (s) and G m (s) have the equal number of poles s i with Res i ≥ 0 [39,50].

Lemma 9. The control system is robustly stable if and only if
where the nominal complementary sensitivity function reads with the nominal parameters in Table 4.

Lemma 10. The control system satisfies inequality (35) and robust stability condition (42) if and only if
where S 0 is given by (45) Following the design step from the preceding subsection, uncertainty weight functions W M (s) are given by (40) as well since they do not depend on the used control system structure. Contrariwise, the choice of W P (s) can be different. When adopting its form (41), Table 7 is eventually obtained. It is worth noting that m 0 = m 1 = a (see Table 6) has been set again, and the weighting parameter midvalue λ = 0.5 in controller (30) Fig. 15(a-c).
The It is worth noting that the robust performance problem with the Best-fit model cannot be solved by a less conservative selection of W P (s) as the problem is caused by the robust-stability term in (44).
Again  Some other parametric space regions have been computationally explored again, providing better performance values according to (44); however, due to a high control action and a high-speed control action are not suitable for the eventual control experiments. Note that particular robust performance measure values are similar to those for λ = 0.25 and λ = 0.5.
Based on the analysis above, the parameters of controller (30)- (31) provided to the reader in Table 8 have eventually been chosen for further laboratory control experiments.

Laboratory control experiments
Experimental verification of robust controllers designed in Sections 4 and 5 for the 1DoF and TFC control systems follows. The nominal case and several perturbations are considered. The received responses are then evaluated using some performance measures.
Denote controlled plant inputs and outputs by u(t) = Δ s P H (t), y(t) = Δ s ϑ CO (t), respectively, for simplicity. That is, the zero inputoutput values agree with the steady state (14).

Control responses for 1DoF
The nominal responses of u(t) and y(t) are displayed in Fig. 19(a and b). Now, consider a perturbation of u P (t) = 4 V (instead of the nominal value u P (t) = 5 V). Then, the corresponding control responses are provided in Fig. 20(a and b). Note that such a perturbation significantly impacts process delays. Let another process perturbation be u C (t) = 4 V (instead of the nominal value u C (t) = 3 V), see Fig. 21(a and b). Finally, let us test a decreased ambient temperature ϑ a (t) = 21 • C (instead of the nominal value ϑ a (t) = 24 • C) ensured by a room thermostat; see Fig. 22(a and b).

Control responses for TFC
The nominal responses of u(t) and y(t) for λ = 0.3 and λ = 0.7 are displayed in Fig. 23(a and b) and Fig. 24(a and b), respectively. Responses for the perturbed case u P (t) = 4 V with λ = 0.3 and λ = 0.7 are given in Fig. 25(a and b) and Fig. 26(a and b), respectively. Control performance under perturbation u C (t) = 4 with λ = 0.3 and λ = 0.7 are displayed in Fig. 27(a and b) and Fig. 28(a and b), respectively. Finally, control responses under the ambient temperature perturbation ϑ a (t) = 21 • C for λ = 0.3 and λ = 0.7 are displayed in Fig. 29(a and b) and Fig. 30(a and b), respectively.

Experimental results evaluation
Several metrics are used for the evaluation of the control responses. Besides the IAE (12), output temperature overshoots/undershoots (Δy max ) are measured, and three forms of the total variation (TV) are computed [55,56] to evaluate the monotonicity and deviation from it.
Controlled outputs are subject to an adjusted TV (denoted by TV 1 here in (46)) that expresses the difference between the total signal path and the minimum possible path required to change the signal from the initial state y k0 to the final state y k1 : It holds that the value of TV 1 (y) equals zero for monotonic outputs.
The basic TV (denoted by TV 0 here in (47)) is used to control action u(t) in those intervals, where the reference signal is linear-wise, as another one is applicable for a step-wise r(t), where u extr is an extreme value lying u extr ∕ ∈ (u k0 , u k1 ), i.e., outside the interval given by the initial and final values. TV 2 (u) in (48) expresses the deviation from an ideal single-pulse waveform (composed of two monotonic intervals) based on the Feldbaum theorem [55]. TVs for the control action mean the control effort, which, i.a., corresponds to the actuators' lifetime.
In addition, the integral of P H (t) is evaluated since (49) represents energy consumption, which is closely related to a very topical issue of sustainability and energy saving.
To evaluate the effects of step-wise and linear-wise reference changes, their increase and decrease, and the impact of load disturbance separately, let us calculate the performance measures for five disjoint time intervals: Then, particular k 0 and k 1 correspond to the initial and final values, respectively, of intervals (50). Notice that I 5 serves only to the evaluation of the disturbance effect. Besides, the total IAE and EC for y(t) and P H (t), respectively, are calculated.
The results are summarized in Table F1 (1DoF), Table F2 (TFC, λ = 0.3), and Table F3 (TFC, λ = 0.7) that can be found in Appendix. Note that the values of EC(P H ) are in kWh in the tables.
Let us select some distinct observations from the data. Regarding the 1DoF control system, relatively high values of m 0 or m 1 (see also Table 6) yield a chattering of u(t), which results in enormous values of TV 0 (u) and TV 2 (u), and also higher overshoots/undershoots. Model 1 and the Best-fit models give low IAEs in the nominal case, but the latter causes a high overshoot after the reference step change. The Best-fit model also provides the best output reaction to the load disturbance (IAE and Δy max ). On the other side, this model is susceptible to perturbations (compared to the other models), which confirms the hypothesis that an almost exact nominal model might not be useful for robust control design and real-world control under uncertainties. The change of u P represents the most distinctive perturbation, mainly due to the change of delays (see Table 3). On the contrary, the responses are almost insensitive to an ambient temperature change. The data shows that Model 1 provides excellent responses under perturbations. Surprisingly, despite high IAEs given by low m 0 , m 1 , the Relay model proves to be sufficient enough for control objectives. An unpleasant effect of perturbations and "fast" controller settings is the existence of the wind-up effect on u(t) that reaches its physical limits. Notice that the Best-fit model causes the wind-up even in the nominal case ( Fig. 19(a,b)). There exist several principles to tackle this problem, e.g., a functional approach for time-delay systems [29]. This problem represents one of the tasks of our future research.
Regarding the TFC control system with λ = 0.3, IAEs for the Original model give better values compared to 1DoF for step-wise reference changes; however, it does not hold for other models. Contrariwise, Model 1 and the Relay model yield lower IAEs when linear-wise reference tracking. A significant undesirable effect of TFC is an increase in overshoots/undershoots after a step change of r(t). Hence, the suggestion here is to use a filtered reference signal without abrupt changes in practice. In the cases of the Original model and Model 1, a better response to the load disturbance can be observed. An interesting observation can be made on the energy consumption: Its value seems almost invariant to the used model and control system and depends solely on a particular perturbation.
The setting λ = 0.7 brings about higher IAEs but lower overshoots after the reference change compared to λ = 0.3. Contrariwise, IAEs given be a reaction to the disturbance are better. Overall, the TFC control system cannot be indicated as better compared to the 1DoF structure. However, the possibility of tuning λ value brings more options when finding a trade-off between reference tracking and disturbance rejection.

Link to the Smith predictor
Finally, let us concisely compare to the well-establish Smith dead-time compensator [51]. As the process model structure (16)- (17) includes the input-output delay, such a comparison can be made. We also point out a significant drawback of this approach.
The Smith predictor structure is depicted in Fig. 31.
In the scheme, G m0 (s) means a process model transfer function without the input-output delay. This modelled delay (i.e., its estimation) is represented by block e − τms . The Smith-predictor controller has transfer function C Sm (s). Two research questions arise: 1) What is the relation between C Sm (s) and C R (s), C Q (s)? 2) Does C Sm (s) work well?
For the 1DoF control system, the following relation can be derived based on the matching of reference-to-output transfer functions (i.e., T(s)) See Appendix F for a sketch of the proof. Note that controller C R (s) is given by (23)- (24) For the TFC control system, the analogous matching yields The reader is referred to Appendix F again. Transfer functions of controllers C R (s), C Q (s) for the heating-cooling process in question  can be found in (30)- (31). It is worth noting that whereas (51) is independent of the actual process dynamics, controller (52) also depends on G(s). It implies that the agreement of TFC and the Smith predictor depends on an accurate estimate of the process dynamics.
Settings (51) and (52), however, bring about a significant drawback of this design: Theorem 1. The use of the Smith predictor with controllers (51) and (52) for process models (23)- (24) and (30)- (31), respectively, cannot guarantee an asymptotic linear-wise load disturbance rejection. Proof of Theorem 1 can also be found in Appendix G. It is worth noting that the designed 1DoF control system rejects the linear-wise disturbance; however, TFC does not, as indicated in the appendix as well.

Conclusions
A detailed modelling, identification, and control study has been presented in this paper. A laboratory looped heating-cooling system with significant input-output and internal delays has been considered as the process in question. A thorough revision and reformulation of the model parameters identification procedure have been designed, resulting in a novel process model reflecting process and measurement uncertainties and perturbations. Besides, an accurate nominal model has been assembled as well for comparison. The control design has included a controller structure derivation for two different control systems, namely the 1DoF and the TFC structures, and the application of robust stability and robust performance conditions in detail. As a result, suitable controller parameter settings have been obtained. Laboratory experiments have proven the applicability of the controllers proposed based on the   robust model in the nominal case and for several process perturbations. Moreover, comparisons to some other alternative models have been made. Namely, a controller designed based on a model obtained from a recent relay-feedback experiment has also proven to be capable of controlling the process under perturbations and a load disturbance. Contrariwise, a controller that best matches the nominal process data could not meet the robust performance conditions. A link to the well-established Smith dead-time compensator has concisely been introduced, and issues with process model estimation and disturbance attenuation have been highlighted.
The main findings of this research can be summarized as follows:   1) The static-parameter approximation error of the proposed model (Model 1) is less than the Original model.
2) Operating ranges of process delays have been estimated based on the dynamic-parameter identification for all considered models.
3) The Relay model has evinced the worst dynamic responses when matching the measured data. 4) The Best-fit model could not meet the robust stability and robust performance conditions for both control system structures.
Model 1 has been very close to the borders for the 1DoF control system.   5) The Original and Relay models have satisfied the robust performance condition for slow control response settings in the TFC structure. 6) The Original model has given better step-wise reference with the TFC structure than the 1DoF scheme. 7) The Best-fit model matching the nominal step responses most accurate has provided poor robust control performance and could not be used for real-world control under uncertainties. 8) Relay model has proven to be sufficient enough for control purposes. 9) Model 1 has resulted in a lower IAE when linear-wise reference tracking with the TFC control system. 10) TFC has improved a load disturbance response but increased overshoots/undershoots after a reference step change compared to 1DoF. 11) The change of the pump input voltage has represented the most distinctive perturbation, mainly due to the change of delays. 12) Overall energy consumption has remained almost invariant to the used model and control system and depended only on perturbations. 13) The value of λ can serve as a tuning knob to reach a trade-off between reference tracking and disturbance rejection. 14) The analogous Smith predictor could not satisfy asymptotic linear-wise load disturbance rejection.
Possible future research can concern overshoot reduction (e.g., based on a model predictive control framework), anti-wind-up protection, a multivariable control design [57], or the use of a multiloop control system (with the auxiliary controlled variable), or the implementation of advanced parameter optimization techniques when identification [58]. In any case, we are motivated by a constant effort to create advanced modelling and control techniques to reduce energy dependence while maintaining sufficient user comfort.

Author contribution statement
Libor Pekař: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.
Radek Matušů: Conceived and designed the experiments; Analyzed and interpreted the data; Wrote the paper. Petr Dostálek: Performed the experiments; Wrote the paper. Mengjie Song: Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Data availability statement
Data will be made available on request.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Libor Pekar reports financial support was provided by Tomas Bata University in Zlin. Libor Pekar reports financial support was provided by College of Polytechnics Jihlava. Libor Pekar reports a relationship with Tomas Bata University in Zlin that includes: employment.

Appendix B
Consider the difference output representation (8) and for the operation point (14), and introduce the notation Let ϑ a (t) = ϑ a = s ϑ a = 24 = const., which i.a. means that the ambient temperature is considered to be a model parameter with the given nominal value, not a disturbance input. This option is further used when performing robustness tests.
Apply the simple linearization rule for DDEs (1), (3), and (4) after substituting the algebraic formulae (2), (5), (6), and the delayed relation (7), where Δϑ j (t) stands for entries of Δϑ(t) and Δu i (t) that of Δu(t). Symbol Ω represents the operating point (14). Then, the following linearized state-space model is received [26] d dt in which the matrices read where thwith static gains adjustmente entries are , ] , where The reader is referred to Ref. [26] for detailed forms of particular transfer functions in G(s).

Appendix C
A concise description of the R QM ring follows [49,50]. where s ∈ C is the Laplace variable, q ij ∈ R are coefficients, 0 = τ i0 < τ i1 < … < τ iki ∈ R mean delays, and q n (s) is the so-called associated exponential polynomial. Denote by r Q [s] the set of quasi-polynomials.
Definition C.2 A fraction F(s) = q num (s)/q den (s) where q num (s), q den (s) ∈ R Q [s] is formally stable if q n (s) in q den (s) has no zero in the left half of the s-plane. |F(s)| < ∞.

Appendix D
Consider the process submodel transfer function (17)   the numerator of which is given by (24). According to (20) software  TFC  two-feedback-controllers  TTL  transistor-transistor logic  TV  total variation  1DoF one-degree-of-freedom